Motivation: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

• For ease of exposition, dropping education.

Motivation: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

• Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]

• Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]

Motivation: Mincer Wage Equation

• Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

• Let \(\overline{\mbox{exp}}_N\) denote sample mean of experience.

• Estimate average effect of experience by \[\widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N,\]

  • Estimator is a non-linear function of \((\widehat \beta_{1,N}, \widehat \beta_{2,N}, \overline{\mbox{exp}}_N)\).

Motivation: Mincer Wage Equation

• Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

• Estimate experience level that maximizes log wage by \[\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right|.\]

  • Estimator is a non-linear functions of \((\widehat \beta_{1,N}, \widehat \beta_{2,N})\).

Implement Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))

• Consider estimation using CPS 1985 data.

  • available in package AER.

  • after loading AER, see ?CPS1985 for description of data.

• Consider only using observations with \(12\) years of schooling.

Implement Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))

[1] "Mean Experience: 18.14"

Call:
lm(formula = I(log(wage)) ~ experience + I(experience^2), data = df)

Coefficients:
    (Intercept)       experience  I(experience^2)  
      1.5829802        0.0352680       -0.0005831  

Implement Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp

[1] "Estimated Avg Effect Experience: 0.0141"

[1] "Estimated Experience that Maximizes Wage: 30.2"

Implement Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp

Motivation: Mincer Wage Equation

• We have not covered methods for constructing standard errors or performing inference in this context.

• How to compute standard errors? construct confidence intervals? perform hypothesis tests?

Traditional approach: Delta Method.

• Asymptotic approximation based on linear approximation from Mean Value Theorem combined with LLN, CLT.

• Limitations:

  • only justified asymptotically,
  • sometimes provides poor approximation,
  • derivations using delta method can be tedious, chance of algebraic errors.

• we will not cover in Econ 123

Use Monte Carlo Simulation?

• If knew distribution of population being sampled, could use Monte Carlo simulation drawing from that distribution to find distribution of our estimator.

  • But we don’t know population distribution.

Bootstrap principle

• Treat sample of size \(N\) as if it were the population.

• Perform MC simulation, creating many bootstrap samples of size \(N\) by resampling with replacement from original sample.

• Compute estimator on each sample, use distribution of estimates across samples to approximate distribution of estimator.

• Difference between drawing from population vs from sample should vanish as \(N \rightarrow \infty\).

Bootstrap: Basic Principle

Bootstrap Method

• Limitations:

  • only justified asymptotically,
  • sometimes provide poor approximation,
  • can be computationally intensive

• Advantages:

  • Can bypass derivations, complicated formulas. Less tedious to implement. Removes chance of algebraic error in derivations.
  • Sometimes provides better approximation than delta method.

Bootstrap Method

• There are many versions of the bootstrap, but we will focus on nonparametric bootstrap and corresponding bootstrap:

  • estimate of bias,
  • standard errors,
  • normal-approximation CI,
  • percentile CI,
  • BCa: bootstrap accelerated bias-corrected percentile CI.

Sample and Estimator

• Suppose \(\mathbf{W} = (W_1,W_2,...,W_N)\) is observed sample
  • e.g., \(W_i = (Y_i,X_i)\) in regression context.
• Estimator \(\widehat \theta = s(\mathbf{W})\)
  • expressing estimator as function of data, rule for how to use data to compute estimate.
  • in our application, \(\widehat \theta\) is estimator of average effect of experience or of experience level that maximizes wages.

Example: Bootstrap Estimate Bias

• Given \(B\) bootstrap replications, define the boostrap estimate of bias by \[\widehat{\mbox{bias}}^*= \frac{1}{B} \sum_{b=1}^B \widehat \theta^*(b) - \widehat \theta \] where \(\widehat \theta\) is estimate from original sample.

• why centered at \(\widehat \theta\)?

Example: Algorithm for Percentile CI

• Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\)

  • compute \(q^*_{\alpha}\) as the value such that fraction \(\alpha\) of the bootstrap replications are less than \(q^*_{\alpha}\).

• Bootstrap Percentile CI: \[C^{PC}= [q_{\alpha/2},~ q_{1-\alpha/2}].\]

Pros and Cons Bootstrap Methods

• Bootstrap standard errors:
  • Pros:
    • intuitive, easy to implement,
    • researchers familiar with s.e., is the norm to report s.e..
  • Cons:
    • Valid only under relatively strong conditions,
    • May be erratic, poorly behaved.

Pros and Cons Bootstrap Methods

• Normal-approximation bootstrap CI:

  • Pros:
    • Intuitive, easy to implement,
    • Researchers familiar with normal CIs,
    • Many estimators are asymptotically normal.

Pros and Cons Bootstrap Methods

• Normal-approximation bootstrap CI:

  • Cons:
    • Invalid if bootstrap s.e. are invalid,
    • Poorly behaved if bootstrap s.e. are poorly behaved,
    • Normality may be bad approximation
      • especially poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

• bootstrap percentile CI:

  • Pros:

    • intuitive, easy to implement,

    • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

  • Cons:

    • can provide poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

• BCa:

  • “bootstrap accelerated bias-corrected percentile interval”

  • adjusts percentile CI for median bias and skewness in bootstrap distribution

Pros and Cons Bootstrap Methods

• BCa:

  • Pros:
    • provides more accurate CIs than above options, especially when estimator is biased or has skewed distribution.
  • Cons:
    • complicated, not intuitive.
    • (somewhat) more computationally intensive.

Implement: Avg Effect of Experience

• We now implement bootstrap s.e. and CIs for average effect of experience.

  • \(\theta = \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}[\mbox{exp}],\)

  • \(\widehat \theta = \widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N.\)

• We first proceed in base R, then consider using package boot.

• First step: create bootstrap samples.

Implement: Avg Effect of Experience

• Create bootstrap sample:

  • Randomly draw \(N\) integers, with replacement, from \((1,...,N)\), with each integer equally likely to be draw,

  • Create new sample formed by observations with those indices.

  • Some observations from original sample will be drawn multiple times into bootstrap sample, others will not appear in bootstrap sample.

N <- nrow(df) 
indices.b1<-sample(1:N,replace=TRUE)  # draw indices with replacement from (1,...,N), N sample size
tab.b1 <- df[indices.b1,] # construct bootstrap sample

Bootstrap Replications Avg. Effect Experience

• given bootstrap sample, \(\mathbf{W}^{*b}=(W^{*b}_1,W^{*b}_2,…,W^{*b}_N)\), compute \(\widehat \theta^*(b)=s(\mathbf{W}^{*b}).\)

  • compute mean experience on bootstrap sample, \(\overline{\mbox{exp}}^{~*b}\).
  • run regression on bootstrap sample, obtain bootstrap sample OLS coefficients \((\widehat \beta_{0}^{*b},\widehat \beta_{1}^{*b}, \widehat \beta_{2}^{*b})\)
  • compute \[\widehat \theta^*(b) = \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b}.\]

Bootstrap Replications Avg. Effect Experience

#  tab.b1 <- df[indices.b1,], bootstrap sample
mean.exp.b1 <- mean(tab.b1$experience) #  sample mean experience on bootstrap sample
reg.b1 <- lm(I(log(wage))~ experience+I(experience^2),data=tab.b1) # OLS on bootstrap sample
# form bootstrap replication of estimated mean effect experience
a <- c(0,1,2*mean.exp.b1)
theta.b1 <- as.numeric(t(a)%*%coef(reg.b1)) 

Construct \(B\) Bootstrap Replications

• We have created 4 bootstrap samples, and corresponding 4 bootstrap replications.

• We want to create \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications, for \(B\) large, for example, \(B=10,000\).

• How to repeat steps many times?

  • we first consider approach in base R using replicate,
  • we then consider using functions from package boot.

Use replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• Create user defined function and use replicate to generate \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• User defined function takes two arguments:
  1. original sample as data frame,
  2. indices;
• Function returns bootstrap replication.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))

• Inputting user defined function, original data and random sample of indices into replicate function, which is returning \(B\) bootstrap replications.

• Use output to compute bootstrap s.e. or CI.

BCa Confidence Interval

Bootstrap accelerated bias-corrected percentile interval (BCa):

• improves upon bootstrap percentile CI by adjusting for bias and skewness in sampling distribution.

• method not intuitive, is complicated, justification is more advanced.

• easy to implement using boot.ci.

• should use it when there is substantial bias/skewness in sampling distribution.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• user defined function takes two inputs:
  • data frame (original sample)
  • indices;
• function returns bootstrap replication.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• boot taking as inputs:

  • original sample (data frame),
  • user defined function,
  • number of replications,

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot1, R = B, stype = "i")


Bootstrap Statistics :
      original        bias    std. error
t1* 0.01410963 0.00004975603 0.002904684

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

• boot.ci taking as inputs:

  • saved output from boot,
  • choice of confidence level,
  • type of bootstrap CI.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0084,  0.0198 )   ( 0.0085,  0.0199 )   ( 0.0083,  0.0197 )  
Calculations and Intervals on Original Scale

Rerun Analysis?

• Can rerun bootstrap simulation to see if results are stable across simulations.
• If not?
  • insufficient number of replications? increase \(B\)?
  • problem with bootstrap method in particular example?
    • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.

Rerun Analysis?

• Can rerun bootstrap simulation to see if results are stable across simulations.
• If not?
  • problem with bootstrap method in particular example?
    • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.
    • percentile and BCa CIs are valid under weaker assumptions and tend to be more robust/stable.

Implement: Exp Level Max. Wages

• We now implement bootstrap s.e. and CIs for experience level that maximizes log wage.

  • \(\theta = \left | \frac{\beta_{1}}{2 \cdot \beta_2} \right|,\)
  • \(\widehat \theta = \left | \frac{\widehat \beta_{1}}{2 \cdot \widehat \beta_2} \right|,\)

• Key difference from before:

  • parameter involves ratio of coefficients,
  • estimator divides by estimated coefficient,
  • cause for worry?

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot2, R = B, stype = "i")


Bootstrap Statistics :
    original   bias    std. error
t1* 30.23942 4.237817    142.7423

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-2013.71,  2042.13 )   (   25.29,    51.45 )   (   25.29,    51.23 )  
Calculations and Intervals on Original Scale

Paired vs Residual vs Wild Bootstrap

• We have done nonparametric bootstrap for OLS, also called “pairs bootstrap” in regression context.

• Other bootstrap procedures for regressions:

  • Residual bootstrap,

  • Wild bootstrap.

• Tradeoffs?

Non-i.i.d data

• We have done bootstrap for i.i.d. data.

• Can modify bootstrap procedure for

  • clustered sampling
    • straight forward modification, resample clusters
  • time series
    • less straight forward (block bootstrap).

  .